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Abstract. We study a spatially explicit harvesting model in periodic or bounded environ- 
ments. The model is governed by a parabolic equation with a spatially dependent nonlinearity 
of Kolmogorov-Petrovsky-Piskunov type, and a negative external forcing term —8. Using sub- 
and supersolution methods and the characterization of the first eigenvalue of some linear ellip- 
tic operators, we obtain existence and nonexistence results as well as results on the number of 
stationary solutions. We also characterize the asymptotic behavior of the evolution equation as 
a function of the forcing term amplitude. 

In particular, we define two critical values 8* and 82 such that, if 8 is smaller than S* , 
the population density converges to a "significant" state, which is everywhere above a certain 
small threshold, whereas if 8 is larger than 82, the population density converges to a "remnant" 
state, everywhere below this small threshold. Our results are shown to be useful for studying 
the relationships between environmental fragmentation and maximum sustainable yield from 
populations. We present numerical results in the case of stochastic environments. 

1. Introduction 

Overexploitation has led to the extinction of many species [3]. Traditionally, models of or- 
dinary differential equations (ODEs) or difference equations have been used to estimate the 
maximum sustainable yields from populations and to perform quantitative analysis of harvest- 
ing policies and management strategies [17j . Ignoring age or stage structures as well as delay 
mechanisms, which will not be treated by the present paper, the ODEs models are generally of 
the type 

(1.1) ^ = F{U )-Y(U), 

where U is the population biomass at time t, F(U) is the growth function, and Y(U) corresponds 
to the harvest function. In these models, the most commonly used growth function is logistic, 
with F(U) = U{n — vU) (see [5], [25], |35|). where \i > is the intrinsic growth rate of the 
population and v > models its susceptibility to crowding effects. 

Different harvesting strategies Y(U) have been considered in the literature and are used in 
practical resource management. A very common one is the constant-yield harvesting strategy, 
where a constant number of individuals are removed per unit of time: Y(U) = 5, with 5 a 
positive constant. This harvesting function naturally appears when a quota is set on the har- 
vesters [31] , [32] , [38] . Another frequently used harvesting strategy is the proportional harvesting 
strategy (also called constant- effort harvesting), where a constant proportion of the population 
is removed. It leads to a harvesting function of the type Y(U) = 5U. 
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Much less has been done in this field using reaction-diffusion models (but see [23], [26], |29j). 
The aim of this paper is to perform an analysis of some harvesting models, within the framework 
of reaction-diffusion equations. 

One of the most celebrated reaction-diffusion models was introduced by Fisher [15] and Kol- 
mogorov, Petrovsky, and Piskunov [22] in 1937 (we call it the Fisher-KPP model). Since then, it 
has been widely used to model spatial propagation or spreading of biological species into homo- 
geneous environments (see books [25], [28], and [30j for a review). The corresponding equation 
is 



where u = u(t, x) is the population density at time t and space position x, D is the diffusion 
coefficient, and [i and v still correspond to the constant intrinsic growth rate and susceptibility 
to crowding effects. In the 1980s, this model was extended to heterogeneous environments by 
Shigesada, Kawasaki, and Teramoto |37J. The corresponding model (which we call the SKT 
model in this paper) is of the type 

(1.3) Ut = DX7 2 u + u(fi(x) — u{x)u). 

The coefficients n{x) and v{x) now depend on the space variable x and can therefore include some 
effects of environmental heterogeneity. More recently, this model revealed that the heterogeneous 
character of the environment plays an essential role in species persistence, in the sense that for 
different spatial configurations of the environment a population can survive or become extinct, 
depending on the habitat spatial structure [8], [12], [34], [36] . 

As mentioned above, the combination of a harvesting model with a Fisher-KPP population 
dynamics model, leading to an equation of the form u% = DX7 2 u + it(/x — vu) — Y(x,u), has 
been considered in recent papers, either using a spatially dependent proportional harvesting 
term Y(x,u) = q{x)u in [26 , [29J, or a spatially dependent and time-constant harvesting term 
Y(x) = h{x) in [23]. In these papers, the models were considered in bounded domains with 
Dirichlet (lethal) boundary conditions. 

Here we study a population dynamics model of the SKT type, with a spatially dependent 
harvesting term Y(x,u): 



We mainly focus on a "quasi-constant-yield" case, where the harvesting term depends on u only 
for very low population densities (ensuring the nonnegativity of u). We consider two types of 
domains and boundary conditions. In the first case, the domain is bounded with Neumann 
(reflective) boundary conditions; this framework is often more realistic for modeling species 
that cannot cross the domain boundary. In the second case, we consider the model (jl.4p in the 
whole space M. N with periodic coefficients. This last situation, though technically more complex, 
is useful, for instance, for studying spreading phenomena [7], [9], and for studying the effects 
of environmental fragmentation, independently of the boundary effects. Lastly, note that the 
effects of variability in time of the harvesting function will be investigated in a forthcoming 
publication [13] . 

In section 2, we define a quasi-constant-yield harvesting reaction-diffusion model. We prove, 
on a firm mathematical basis, existence and nonexistence results for the equilibrium equations, as 
well as results on the number of possible stationary states. We also characterize the asymptotic 
behavior of the solutions of (11.41) . In section 3, we illustrate the practical usefulness of the results 
of section 2, by studying the effects of the amplitude of the harvesting term on the population 
density in terms of environmental fragmentation. Lastly, in section 4, we give new results for 
the proportional harvesting case Y(x, u) = q(x)u. 



(1.2) 



u t = DV 2 u + u(/i — uu) 



(1.4) 



ut = DV 2 u + u([a(x) — u(x)u) — Y(x, u). 
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2. Mathematical analysis of a quasi-constant- yield harvesting 

reaction-diffusion model 

For the sake of readability, the proofs of the results of section 2 are postponed to section 2.5. 

2.1. Formulation of the model. In this paper, we consider the model 

(2.1) u t = DV 2 u + u(p(x) - u{x)u) - 5h(x)p £ (u), (t, x) G M. + x Q. 

The function u = u(t, x) denotes the population density at time t and space position x. The 
coefficient D, assumed to be positive, denotes the diffusion coefficient. The functions p{x) and 
v{x) respectively stand for the spatially dependent intrinsic growth rate of the population, and 
for its susceptibility to crowding effects. Two different types of domains Vt are considered: either 
Q = M. N or Q is a smooth bounded and connected domain of M. N (N > 1). We qualify the first 
case as the periodic case, and the second one as the bounded case. In the periodic case, we 
assume that the functions p(x), v{x), and h(x) depend on the space variables in a periodic 
fashion. For that, let L = (L\, . . . ,Ljf) G (0, +00)^. We recall the following definition. 

Definition 2.1. A function g is said to be L-periodic if g{x+k) = g(x) for all x = (x\, . . . ,xn) £ 
R N and k G L{L x ••• x L N Z. 

Thus, in the periodic case, we assume that p,, u, and h are L-periodic. In the bounded case 
we assume that Neumann boundary conditions hold: lp = on d£l, where n is the outward 
unit normal to dQ. The period cell C is defined by 

C:= (0,Li) x ••• x (0,L N ) 

in the periodic case, and in the bounded case we set 

c -.= n, 

for the sake of simplicity of some forthcoming statements. 
We furthermore assume that the functions /U and v satisfy 

(2.2) /i,i/GL°°(0) and 3 v ,V G R s.t. < v < v(x) < V V x G ft. 

Regions with higher values of fJ-(x) and lower values of v(x) will be qualified as being more 
favorable, while, on the other hand, regions with lower (J,(x) and higher v{x) values will be 
considered as being less favorable or, equivalently, more hostile. 

The last term in (|2.ip . 5h{x)p £ {u), corresponds to a quasi-constant-yield harvesting term. 
Indeed, the function p £ satisfies 

(2.3) p e G C 1 (M), p' £ > 0, p £ (s) = Vs < and p £ (s) = 1 Vs > e, 

where e is a nonnegative parameter. With such a harvesting function, the yield is constant in 
time whenever u > e, while it depends on the population density when u < e. In what follows, 
the parameter e is taken to be very small. As we prove in the next sections, there are many 
situations where the solutions of the model always remain larger than e. For these reasons, we 
qualify our model as a quasi- constant-yield harvesting SKT model, the "dominant" regime being 
the constant-yield one. Note that the function p £ ensures the nonnegativity of the solutions of 
()2.ip . From a biological point of view, e can correspond to a threshold below which harvesting is 
progressively abandoned. Considering constant-yield harvesting functions without this threshold 
value would be unrealistic since it would lead to harvest on zero-populations. 
Finally, we specify that 5 > and that h is a function in L°°(£l) such that 

(2.4) 3 a > with a < h(x) < 1 Vx G O. 

We call h the harvesting scalar field, and 5 designates in this way the amplitude of this field. 
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Before starting our analysis of this model, we consider the no-harvesting case, i.e., when 5 = 0. 
We recall the main known results in this case. These results will indeed be necessary for the 
analysis of the quasi-constant-yield harvesting SKT model. 

2.2. The no- harvesting case. When 5 = in (|2.ip . our model reduces to the SKT model 
described by (jl.3p . The behavior of the solutions of this model has been extensively studied in 
[8] and [9]. 

Results are formulated in terms of first (smallest) eigenvalue Ai of the Schrodinger operator 
Cfj, defined by 

£ M := -DV 2 - fi(x)I, 

with either periodic boundary conditions (on the period cell C) in the periodic case or Neumann 
boundary conditions in the bounded case. This operator is the linearized one of the full model 
around the trivial solution. Recall that Ai is defined as the unique real number such that there 
exists a function <p > 0, the first eigenfunction, which satisfies 

, v / -DV 2 (/) - n(x)</> = \i<f> in C, 

1 ' \ <P > in C, |H|oc = 1, 

with either periodic or Neumann boundary conditions, depending on f2. The function <j> is 
uniquely defined by (I2.5j) [7J and belongs to W 2 ' T {C) for all 1 < r < oo (see [I] and [2] for 
further details). We set 

<f) := mm(j)(x). 

We recall that a stationary state p of (|1.3p satisfies the equation 
(2.6) - DV 2 p = p(ji(x) - v{x)p) . 

The following result on the stationary states of (12.6|) is proved in [8]. 

Theorem 2.1. (i) If \\ < 0, then (|2.6j) admits a unique nonnegative, nontrivial, and bounded 
solution, pq. 

(ii) If Ai > 0, the only nonnegative and bounded solution of (|2.6p is 0. 

Moreover, in the periodic case, the solution po is L-periodic. Throughout this paper, po always 
denotes the stationary solution given by Theorem 12. H i. 

In order to emphasize that this solution can be "far" from (see Definition 12.21 and the 
commentary following (|2.10p ). we give a lower bound for p$. 

Proposition 2.1. Assume that \\ < 0; then p$ > — in f2. 

The asymptotic behavior of the solutions of (|1.3p is also detailed in [8j. It is proved that Ai < 
is a necessary and sufficient condition for species persistence, whatever the initial population u° 
is, as follows. 

Theorem 2.2. Let u° be an arbitrary bounded and continuous function in Q such that u° > 0, 
u° ^ 0. Let u(t,x) be the solution of (jl.3p . with initial datum u(0,x) = u°(x). 

(i) If Ai < 0, then u(t, x) — > po(x) in W,'"[ (O) for all 1 < r < oo as t — > +oo (uniformly in 
the bounded case). 

(ii) // Ai > 0, then u(t, x) — > uniformly in £1 as t — > +oo. 

The situation (i) corresponds to persistence, while in the case (ii) the population tends to 
extinction. In what follows, unless otherwise specified, we therefore always assume that Ai < 0, 
so that the population survives, at least when there is no harvesting. We are now in position 
to start our main analysis of steady states and related asymptotic behavior of the solutions of 

(EU). 
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Figure 1. The logistic growth function U t— >■ U{p — vU) (solid line), and the har- 
vesting function U t— >■ 5p £ (U) for three values of 5 (dashed lines). The abscissae of 
the points of intersection of the solid and dashed lines correspond, respectively, 
to remnant (if smaller than e) and significant (if strictly larger than e) station- 
ary solutions of (|2.9p . We observe that the number of significant solutions is as 
follows: one if S < k(e) (case 5 = 5 a ); two if k(e) < 5 < p, 2 /(4u) (case 5 = 5b); 
one if 6 = /i 2 /(4i^) (case 6 = 5*); zero if 5 > p 2 /(Av). The number of nonzero 
remnant solutions is zero or more if 8 < k(e) (depending on the shape of p e ); one 
or more if 5 > k(e), since, from (|2.3p . p' e (0) = 0. We assumed here that Eq = e. 



2.3. Stationary states analysis. As is classically demonstrated in finite dimensional dynam- 
ical systems theory and many problems in the infinite dimensional setting (see, e.g., [39 J), the 
asymptotic behavior of the solutions of f)2. 1 1) is governed in part by the steady states and their 
relative stability properties. In that respect, we study in this section the positive stationary 
solutions of f)2. 1 [) . namely the solutions of 

(2.7) - DV 2 p s = ps(p(x) - v(x)ps) - 5h(x)p e (ps), x £ Q, 

in the periodic and bounded cases. When needed, we may write (12. 7} 5) instead of (12.71) . 
Note that, provided ps > e in f2, p§ is equivalently a solution of the simpler equation 

(2.8) - DX7 2 ps = ps{p{x) - u(x)ps) - 5h(x), x G tl. 

This last equation has been analyzed in the case of Dirichlet boundary conditions in [29], in the 
particular case of constant coefficients p and v. 

Because of the type of harvesting function considered here, we are led to introduce the fol- 
lowing definition. 

Definition 2.2. Set sq := ^ > e. We say that a nonnegative function a is remnant whenever 
maxc<7 < £q> whereas it is significant if it is a bounded function satisfying mined > Eq. 



Remark 2.1. The concepts of remnant and significant solutions, as well as the harvesting term 
5h(x)p e (u), are not classical. In order to clarify these notions, we present in Figure [T] a short 
graphical study of the nonspatial model 

(2.9) ^ = U(p-vU)-5p £ (U)=:k(U), t G M + , 

with constant coefficients p, v > 0. 
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Since Eq is assumed to be small in our model, the remnant solutions of (|2,7p correspond to 
very low population densities. On the other hand, significant solutions are everywhere above £q. 
In particular, a constant yield is ensured in that case. In contrast to the ODE case, stationary 
solutions which are neither remnant nor significant may exist, as outlined in the next theorems. 
However, as we will see while studying the long-time behavior of the solutions of the model (|2 . 1 1) . 
they are of less importance (see Theorem 12.61 and section 3). The threshold £q is different from 
e in general. We had to define remnant and significant functions using £q for technical reasons 
(see the proof of Theorem I2.51 ii. equation (|2.27| ). Since e is assumed to be very small, it has 
no implication on the biological interpretation of our results. Moreover, most of our results still 
work when eo is replaced by e. 

Let us now start our analysis of (|2.7p . In what follows, we always assume that 

-Ai0 

(2.10) ,o < , 

so that, in particular, from Proposition l2.lt the solution po °f (|2.6p is significant. 

We begin by proving that there exists a threshold 5* such that, if the amplitude 5 is below 
5* , (I2.7P admits significant solutions, while it does not in the other case. 



Theorem 2.3. Assume that Ai < 0; then there exists 5* > such that 

(i) if 8 < 8* , there exists at least a positive significant solution p$ < pn of ([2? 

(ii) if 5 > 6* , there is no positive significant solution of (|2.7p . 

Remark 2.2. There is no positive bounded solution of (|2,7p whenever Ai > 0. 

Under stronger hypotheses, we are able to prove that (|2.7p admits at most two significant 
solutions. In order to state this result, we need some definitions. Let G be the space defined by 

(2.11) G:=H l (C) 
in the bounded case, and by 

(2.12) G := H^ er = {V G Hl c (R N ) such that ip is L-periodic} 

in the periodic case. Let us define the standard Rayleigh quotient: for all if) G G, ip 0, and 
for all a G L°°(C), 



(2.13) K a ($) :-- 



f c D\Vi/)\ 2 - < j(x)jj 2 



According to the Courant-Fischer theorem (see, e.g., [6]), the second smallest eigenvalue A2 of 
the operator £ M can be characterized by 

(2.14) A2 = min max T^JMi). 

E k CG,dim(E k )=2 4>£E k , ^0 



This characterization is equivalent to the classical one given in [18J . 
We are now in position to state the following theorem. 

Theorem 2.4. Assume that Ai < < A2; then, in the bounded case, (|2.7p admits at most 
two significant solutions. In the periodic case, (j2.7|) admits at most two L' -periodic significant 
solutions for all V G (0, +00)^. Moreover, under these hypotheses, if two solutions pis and p2,s 
exist, they are ordered in the sense that, for instance, pi§ < £>2,<5 $7. 

Remark 2.3. Similar methods also allow us to assess a result on the number of solutions of 
(12. 8j) . Indeed, if Ai < < A2, then we obtain that (|2.8p admits at most two nonnegative bounded 
(and periodic in the periodic case) solutions. If these solutions exist, they are ordered. 
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In the periodic case, Theorem 12,41 also gives some information on the periodicity of the signif- 
icant solutions of (|2.7j) . which are actually found to have the same periodicity as the coefficients 
of (12.71). as seen in the next result. 



Corollary 2.1. Assume that X\ < < A2. Then, in the periodic case, the significant periodic 
solutions of ()2.7p are L-periodic. 

The fact that Ai < is directly related to the instability of the trivial solution in the SKT 
model. The additional condition A2 > in this theorem is linked to the existence of a stable 
manifold or center manifold of the steady state of the SKT model, in some appropriate func- 
tional spaces (see [39]). Therefore, the assumptions of Theorem 12.44 an d the Krein Rutmann 
theory, allow us to conclude that under these assumptions the unstable manifold of is of di- 
mension equal to one or equivalently the stable manifold is of codimension 1. Such results on 
multiplicity of solutions of elliptic nonlinear equations with a source or sink term have been 
investigated in the past and are known nowadays as being of Ambrosetti-problem type. These 
results also involve manifolds of codimension 1 (in the functional space of forcing) and first and 
second eigenvalues (for the Laplace operator only) (see [27] for a survey of these results). 

In any event, Theorem 12.41 relies on the assumption that A2 > 0. In the next proposition, we 
give conditions under which A2 may become positive. 

Proposition 2.2. (i) In the bounded case, if C is a (smooth) domain with diameter d : = 
maxa^ec \\x - y\\ R N, A 2 (C) > D(§ ) 2 - max c /x. 

(ii) In the periodic case, A2(C) > D{j^) 2 — maxc fj-, where Ld denotes the length of the longest 
diagonal of the period cell C. 

For instance, when C = [0, 1] x [0, 1], we have d = Ld = V2; thus, for D = 1 and maxc [i = 4, 
we get A2 > 0.9. However, this lower bound is far from being optimal. Indeed, in all our 
computations of section El and under the same hypothesis on C and D, we always had A2 > 0, 
while maxc fJ. = 10. Sharper lower bounds for A2 can be found in pXJ; however, those bounds 
are also more sensitive to the geometry of the domain and thus less general. They are therefore 
not detailed here. 

We now introduce a result which is important for more applied ecological questions. Indeed, 
one of the main drawbacks of Theorem 12 .31 is that it gives no computable bound for 5* . Obtaining 
information on the value of 5* is precious for ecological questions such as the study of the 
relationships between 5* and the environmental heterogeneities. The next theorem states some 
computable estimates of 5*. 

Let us define 

A?0 A? 
(2.15) 5t := — *~ and 5 2 ~ 1 



Note that neither 5\ nor 62 depend on 5 and e. 

Theorem 2.5. (i) // Ai < and 5 < Si, then there exists a positive significant (and L-periodic 
in the periodic case) solution ps of (|2.7p such that ps > — =^q^y. 

(ii) If Ai < and 5 > 62, the only possible positive bounded solutions of (12.7|) are remnant. 

The lower bound of part (i), for p$, does not depend on e. Thus, there is a clear distinction 
between the remnant and significant solutions. Note that, of course, b\ < 62- 

The formulae fl2. 15f) allow numerical evaluations. An important quantity to compute is the 
size of the gap ^2 — ^1 and its fluctuations in terms of environmental configurations. This question 
is addressed in section [3] through a numerical study. 
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2.4. Asymptotic behavior. In this section, we prove that the quantity 8* in fact corresponds 
to a maximum sustainable yield, in the sense that when 8 is smaller than 8*, the population 
density u(t,x) converges to a significant stationary state of (|2.ip as t — > oo, whereas when 8 is 
larger than 8* , the population density converges to a stationary state which is not significant. In 
fact, when 8 is larger than the quantity 82 defined by (|2,15p we even prove that the population 
converges to a remnant stationary state of (|2.ip . 

We assume here that the harvesting starts on a stabilized population governed by the standard 
SKT model with 8 = 0. From Theorem 12.21 this means that we study the behavior of the 
solutions u(t,x) of our model (|2.ip . starting with the initial datum u(0,x) = po(x). Since we 
have assumed that Ai < 0, it follows from Theorem I2.1| Proposition I2.1| and (|2.10p that pq is 
well defined and significant. 

Let us describe, with the next theorem, the long-time behavior of the population density. 

Theorem 2.6. Let u(t,x) be the solution of f)2. 1 1) with initial datum u(Q,x) = po(x). Then u is 
nonincreasing in t and the following hold: 

(i) If 8 < 8*, u(t,x) — > ps(x) uniformly in £1 as t — >■ +00, where ps is the maximal significant 
solution of (|2.7p . Moreover, p$ is L-periodic in the periodic case. 

(ii) If 8 > 8*, then the function u(i, •) converges uniformly in Q, to a solution of (|2.7p which 
is not significant. 

(iii) If 8 > 82, the function u(t,-) converges uniformly in £1 to a remnant solution of (]2.7p . 

Remark 2.4. If, in addition, we assume that A2 > 0, then Theorem 12.41 says that, whenever 
8 < 8* , (]2.ip admits at most two significant stationary states (which are periodic stationary 
states in the periodic case). In that case, the stationary state ps selected at large times is the 
higher one. If we do not assume that A2 > 0, this stationary state can still be defined as "the 
maximal one" that can be constructed by a sub- and supersolution method (see [3]). 

From the above theorem, we observe that, whenever 8 < 8* , the solution u(t,x) of (12. II) . with 
initial datum po, remains significant for all times t > 0. This ensures a constant yield in time 
and justifies the name of the model. 

Similar results could be obtained for a wider class of initial data. Indeed, with similar methods, 
the convergence of u(t, x) to a significant solution of (12. 7p can be obtained whenever 8 < 8* for 
all bounded and continuous initial data u(0, x) which are larger than the smallest significant 
solution of (12. 7p . In particular, when u(0,x) is larger than the maximal significant solution of 
(12. 7p . u(t, x) converges to this maximal significant solution as t — > +00. A more detailed analysis 
of the basin of attraction related to the maximal significant solution will be further investigated 
in the forthcoming paper [13]. 

Theorem 12.61 shows that the practical determination of 8* is directly linked to the size of the 
gap 82 — 81. As we will see in section [3j this gap (81,82) can be very narrow in certain situations. 
In those cases, the numerical computation of 8\ and 82 therefore gives a sharp localization of 
the maximum sustainable quota 8* E [^1,^2]; which can be of nonnegligible ecological interest. 

2.5. Proofs of the results of section 2. 

Proof of Proposition 12.11 Let <f> be defined by (|2.5p , with the appropriate boundary conditions. 
Set Ko := —= L - Then the function kq4> satisfies 

-DV 2 (K 4>) - n{x)KQ(j) + b>(x)(K (f)) 2 = \lK (f) + h>(x)(K (j>) 2 

= K (f)(\i + v(x)ko4>) < 0. 
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Thus hq4> is a subsolution of (|2,6p satisfied by pq. Since for M G M large enough M is a 
supersolution of (|2.6j) , it follows from the uniqueness of the positive bounded solution po of (I2.6P 
that pq > KQ(f> > — □ 

Before proving Theorem 12.31 we begin with the following lemma. 

Lemma 2.1. For a// 5 > 0, 2/^,5 is a nonnegative bounded solution of (|2.7p . i/ten p,5 < pq. 

Proof of Lemma 12.11 Assume that there exists xq £ Q such that ps(xo) > Po(%o)- The 
function ps satisfies 

-DV 2 ps -ps(n(x) - v{x)ps) = -5h(x)p £ (ps) < 0, 

and thus ps is a subsolution of (|2.6p satisfied by pq. Since for M 6 K large enough M is a 
supersolution of (|2.6p , we can apply a classic iterative method to infer the existence of a solution 
p' of (|2.6p (with Neumann boundary conditions in the bounded case since both pg and M satisfy 
Neumann boundary conditions) such that ps < p' < M. In particular, p' (xo) > po(xo), which 
is in contradiction with the uniqueness of the positive bounded solution of (|2.6p . □ 
Proof of Theorem 12.31 Let us define 

5* := sup{<5 > 0, (|2.7p admits a significant solution}. 

For 5 = 0, we know from Proposition 12.11 that po is a significant solution of (|2.7p . Moreover, 
for 5 large enough, the nonexistence of significant solutions of (|2.7p is a direct consequence of 
the maximum principle (it is also a consequence of the proof of Theorem 12. 51 ii). Thus 5* is well 
defined and bounded. 

Assume that 8* > 0, and let us prove that (|2.71 5*) admits a significant solution. By definition 
of 5* , there exists a sequence (ps k )keN of solutions of (|2.71 5k) with < 5k < 5* and 5k — > 5* as 
k —> +00. Moreover, from Lemma l2.lt e o — Ps k < Po f° r a H k >0. Thus, from standard elliptic 
estimates and Sobolev injections, the sequence (ps k )keN converges (up to the extraction of some 
subsequence) in W^ c , for all 1 < r < 00, to a significant solution ps* of (|2.7| 0"*). 

Now, let < 5 < 5*. Then 

-DV 2 p s * - ps* [n{x) - u(x)p s * ) + 5h(x) = (5 - 5*)h(x) < 0, 

and thus ps* is a subsolution of (|2.71 o"). Since po is a supersolution of (|2.71 5), and < po 5 a 
classical iterative method gives the existence of a significant solution ps of (|2.71 5) (with Neumann 
boundary conditions in the bounded case since both po and ps satisfy Neumann boundary condi- 
tions). This concludes the proof of Theorem 
[2731 □ 

Proof of Theorem 12.41 As a preliminary, we prove that if two solutions exist, then they cannot 
intersect. Let pi } g and p2,<5 be two significant solutions of (|2.7p . In the bounded case, we assume 
that pi : s and p2fi satisfy Neumann boundary conditions. In the periodic case, we assume that 
there exists L' 6 (0, +00)^ such that p\$ and ^2,5 are L'-periodic, and then denote the period 
cell by C . Let us set q$ := P2,s — Pi,8- Then q$ verifies 

(2.16) - DV\ S ~ \n{x) ~ v(x)(pi,s + P2,s)]qs = 0; 
thus, setting p(x) := fi(x) — u(x)(pi i s +P2,s), we obtain 

(2.17) - DV 2 q s - p{x)q s = 0, 

with the same boundary conditions that were satisfied by p\^ and p2,<5- 

Let Ai and A2 be respectively the first and second eigenvalues of the operator C p := —DV 2 —pI. 
Let lZ a {4>), be defined by (|2.13p . Since p(x) < p(x) — 2veq for all x E Q, we get 

Kp{<p) > K^ip) + 2ue 
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for all (p G G", where G' := H l (C) in the bounded case and 

G' := Hp er = {y? G Hi oc (M. N ) such that ip is L'-periodic} 
in the periodic case. Thus, by the classical min-max formula (|2.14p . it follows that 

(2.18) > A 2 + 2v£ Q > 0. 

Furthermore, from (j2.17j) . is an eigenvalue of the operator C p . Thus, (|2.18j) implies that Ai = 0. 
As a consequence, qg is a principal eigenfunction of the operator C p . The principal eigenfunction 
characterization thus implies that qg has a constant sign. Finally, we get that p\ g and p 2y g do 
not intersect each other. 

Let us now prove that (|2.7j) admits at most two significant solutions. Arguing by contradiction, 
we assume that there exist three significant (L'-periodic in the periodic case, for some L' G 
(0, +00)^) solutions pi : g, P2,g, and p^g of f|2.Tj) ■ From the above result, we may assume, without 
loss of generality, that p 3; g > p 2> g > P\.g > e . Set q 2 ,i := p 2 ,s ~ Pi,s and <?3,2 := P3,s ~ P2,s; then 
these functions satisfy the equations 

(2.19) -DV 2 q 2 ,i-p2,i{x)q2,i = 
and 

(2.20) - DV V2 " P3,2 (x)q3, 2 = 0, 

with p 2 ,i := K x ) ~ u ( x )(Pl,5 +P2,s) and p 3t2 := fx(x) - v(x)(p 2t g +P3,<5)- Moreover, g 2 ,l > 
and (/3 5 2 > 0. Thus is the first eigenvalue of the operators C P21 := —DV 2 — p 2 ^I and 
C P3 2 := —DV 2 — P3, 2 I with either Neumann or L'-periodic boundary conditions. 

From the strong maximum principle (see, e.g., [18]) (together with Hopf's lemma in the 
bounded case, and using the L'-periodicity of q 3 ^ 2 in the periodic case), we obtain the existence 
of 6 > such that q 3 ^ 2 > 9. Since the operator C P3 2 is self-adjoint, we have the following formula 

—-3,2 

for its first eigenvalue Ai : 



. — .3 2 

Ai ' = min ftp,, , 2 (</>)• 



Thus 



= ™G'\ Y^P 2 J + - 

— 2,1 
> Ai + 1/0, 

— 2,1 

where Ai is the first eigenvalue of the operator C P21 . Since the first eigenvalues of the operators 
C P2 1 and C P3 2 are both 0, we deduce that > + v9 > 0, hence a contradiction. □ 

Proof of Corollary 12.11 Let pg be a significant L'-periodic solution of (|2.7p . and let G 
ni=i ^iZ. From the L-periodicity of (|2.7p . + k) is also a solution of (|2.7p . By periodicity 
of p,5, the functions and pg(- + fc) intersect each other. Thus, from Theorem 12.41 since pg and 
Ps(' + fc) are both L'-periodic, pg = pg(- + k). Therefore, pg is an L-periodic function. □ 

Proof of Proposition I2.2L In the bounded case, let C be the convex hull of the set C . It was 
proved in [30] that the second Neumann eigenvalue of the Laplace operator — DS7 2 on C was 
larger than D{% ) 2 . Since C C C, we have H X {C) C F 1 (C). Using formula (gjj]) , we thus 
obtain that the second eigenvalue of C p in the bounded case satisfies A2 > D(^) 2 — max^ /J. 
This proves part (i) of Proposition! 
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In the periodic case, since H^ er can be seen as a subset of H l {C), it follows from fl2JH that 

(2.21) A2 > min max TZuNj)- 

E k cm(C)A™{E k )=2i>&E k , ^0 

The period cell C is convex but not smooth enough to assert that the right-hand side of (|2.2ip is 
equal to the second eigenvalue in the bounded case. Let be the longest diagonal of C. Then 
C is included in a ball Bi d of diameter Lj. Thus, from formula (|2.14p . the right-hand side of 
(|2.2ip is larger than the second eigenvalue of on Bi d . From (i), the conclusion of (ii) follows. 
□ 

Proof of Theorem I2.5L part (i). Let Ai and <f> be defined by (|2.5p . and let k be a nonnegative 
real number such that k > £q. Then we have 

-DV 2 (K(j)) - K(/)(n(x) - K(j)v(x)) + 5h(x)p £ (K(j)) < XiK(j) + K 2 (f) 2 u(x) + 5 

(2.22) < K<j>{\i + K(f>v{x)) + 5 

< max{r(Ai + rz/)} + 5, 

where / = {k4>(x), x £ C}. Setting g(j) := r(Ai + Tf), since ||0||oo — 1) and since g is a convex 
function, it follows from (|2.22p that 

(2.23) - DV 2 (K(j)) - K<p(fi(x) - n(f)v{x)) + 8h(x)p E (n(f>) < max{#(K), #(«(/>)} + 5. 

Let us take kq be such that g(no) = gi^ofi), namely ko = ~ ^n^j (note that Kofi > e). We get 

xu 

(2.24) - DV 2 {ko(P) - K^ix) - koM*)) + Sh(x) < j- ~ + 5 < 0, 

from the hypothesis on 5 of Theorem I2.51 i. Therefore, Ko<p is a subsolution of (12.7(1 with either 
L-periodic or Neumann boundary conditions. Moreover, if M is a large enough constant, M is a 
supersolution of (|2.7p with L-periodic or Neumann boundary conditions. Thus, it follows from a 
classical iterative method that there exists a solution ps of (|2.7p . with the required boundary con- 
ditions, and which satisfies kq(J) < ps < M in ft. Moreover, in the periodic case, since Kofi and M 
are L-periodic and since (12. 7j) is also L-periodic, it follows that p# is L-periodic. Theorem I2.51 i is 
proved. □ 

Proof of Theorem 12.51 part (ii). Assume that Ai < 0, 5 > 62, and that there exists a positive 
bounded solution ps of (I2.7P which is not remnant; i.e., 

(2.25) 3 xq with ps(xq) > £q. 

Since fi is bounded from below away from and p$ is bounded, we can define 

(2.26) 7* = inf {7 > 0, jfi > ps in ft} > 0. 

It follows from the definition of 7* that j*fi > ps in ft, and in particular, 7*^(2:0) > Ps(xq) > £q- 
Since ||0||oo = 1, we get 7* > eo- Thus, 

(2.27) 7 *fi > e Q ± = e, 
which implies p £ (-y*<j)) = 1. Thus, h(x) p e (7* fi) > a, and we get 

-L>V 2 (7» - 7* fi(fi(x) - 7*<Ms)) + Sh(x)p £ ^y*fi) > 7*0(Ai + j*<j)v(x)) + 5a 

A 2 

on ft. Moreover, since > and v > u, we have 7*0(Ai + j*fiv(x)) > — Using the fact 
that 5 > ^2 = 7^7, we thus get 

A 2 

(2.28) - Z)V 2 (7*(/ ) ) - 7V(/i(s) - 7*<Ms)) + 8h(x)pe(l*<i>) > 1 + 5a > 

4v 
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on 0. Therefore, is a supersolution of (|2.T|) . Set z := 7*^ — pg. From the definition of 
7*, we know that z > and that there exists a sequence (x n ) n6 N in Q such that z(x n ) — > as 
n —> +00. 

In the bounded case, up to the extraction of some subsequence, i n ->x£Oasn-> +00. By 
continuity, z(x) = 0. Moreover, subtracting (|2.7p from (|2.28p . we get 

(2.29) -DV 2 z + [u{x){^(f) + p s ) + x{x)- n(x)]z > in ft, 

where the function x is defined by x( x ) = ^^( x ) Pe ^\*^)Zp £ J^ X ^ whenever 7*0(x) — p#(x) 7^ 0, 
and x(x) = /^(^(x)) otherwise. Since p s is C , x is bounded. Thus 6(x) := v(x){^*<j) + Ps) + 
x(x) — p{x) is a bounded function. Using the strong elliptic maximum principle, we deduce 
from (|2.29p that z = 0. Thus j*4> = p$ is a positive solution of (|2.7p . It is in contradiction with 
(1^281) . 

In the periodic case, we must also consider the situation where the sequence (x n ) ng pj is not 
bounded. Let (x n ) G C be such that x n — x n G n£i LiL. Up to the extraction of some 
subsequence, we can assume that there exists Xqo G C such that x n — > Xqo as n — > +00. Set 
4>nix) = ^(^ + ^n) and ps,n{x) = ps{x + x n ). From standard elliptic estimates and Sobolev 
injections, it follows that (up to the extraction of some subsequence) ps >n converge in W^'J, for 
all 1 < r < 00, to a function p$ !QO satisfying 

-V 2 (Dp St00 ) - Ps,oc{K x + ^00) - P6,ooV(x + Xqo)) + Sh(x + 5oo)p £ (p5,oo) = 

in M N , while 7*(/>n converges to 7*(/>oo := 7*0(- + ^oo), and 

-V 2 (Z?7*0 OO ) - 7*</>oo(/u(a; + Xqo) - l*(f>oov(x + #oo)) + 5/i(x + Xoo)p e (7*^00) > 

in R . Let us set Zoo(^) := 7*0oo(^) — P<5,oo(^)- Then ^(x) = lim n _ > 4. oc z{x + x n ), and therefore 
Zoo > and Zoo(0) = 0. Moreover, there exists a bounded function such that 

(2.30) - DV 2 z 00 + 600^00 > in R^. 

It then follows from the strong maximum principle that Zoo = 0, and we again obtain a con- 
tradiction. Finally, we necessarily have ps < £0, and the proof of Theorem I2.51 ii is complete. 
□ 

Proof of Theorem 12.61 part (i) . Assume that 8 < 5* . Let p$ be the unique maximal significant 
solution defined in the proof of Theorem 12. 51 i. Then, from Lemma |2.1[ 

(2.31) Ps(x) < po(x) = u(0, x) Vx G fi, 
which implies 

(2.32) p s (x) < u(t,x) in R+ x f2, 

since ps is a stationary solution of ([2.ip . Moreover, since po is a supersolution of ([2.7p . u is 
nonincreasing in time t, and standard parabolic estimates imply that u converges in W?'^ (^), 
for all 1 < r < 00, to a bounded stationary solution of (|2.ip . Furthermore, from ()2.32j) 
we deduce that ps < Uoo < Po- Since p,5 is the maximal positive solution of (|2.7p . it follows 
that tioo = ps- Moreover, in the periodic case, since po and (|2.ip are L-periodic, u(t,x) is also 
L-periodic in x. Therefore the convergence is uniform in £1. Part (i) of Theorem 12.61 is proved. 
□ 

Proof of Theorem 12.61 parts (ii) and (iii). Assume that 5 > 5*. Since is a stationary 
solution of ()2.ip and u(0, x) = po > 0, we obtain that u(t,x) > in R + x Q, and again, from 
standard parabolic estimates, we know that u converges in W, 'J (il) (for all 1 < r < 00) to a 
bounded stationary solution Uqq > of (|2.ip as t — > +00. Moreover, in the periodic case, from 
the L-periodicity of the initial data and of (|2.ip . we know that u(t, •) and «oo are L-periodic. 
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Therefore the convergence is uniform in It follows from Theorem I2.31 ii that cannot be a 
significant solution of (I2.7p . Moreover, if 5 > 62, Theorem I2.51 ii ensures that remnant 
solution of (|2TT) . □ 

3. Numerical investigation of the effects of environmental fragmentation 

We propose here to apply the results of section 2, on the estimation of the maximum sustain- 
able yield, to the study of the effects of environmental fragmentation. A theoretical investigation 
of the relationships between maximum sustainable yield and fragmentation is difficult to achieve 
(see Remark 13. ID . To overcome this difficulty, we propose a numerical study in the case of sto- 
chastic environments. First, we show that the gap 82 — Si, obtained from (|2.15|) and Theorem l2.5l 
remains small whatever the degree of fragmentation is. This gap corresponds to the numerical 
values of the harvesting quota 5 for which we do not know whether the population density will 
converge to a significant or a remnant solution of the stationary equation (12. 7)1 . Second, we show 
that there is a monotone increasing relationship between the maximal sustainable yield S* and 
the habitat aggregation. 

Remark 3.1. In a periodic environment, a simple way of changing the degree of fragmentation 
without changing the relative spatial pattern (favorable area/unfavorable area ratio) is to modify 
the size of the period cell C. Assume that (i(x) = Tj(j^), for some 1-periodic function 77 with 
positive integral and for some L > 0. This means that the environment consists of square cells 
of side L. Setting Xix '■= \\ and 4>l '■= we then have —DA^l — v (-j) <Pl = ^x,L<j>L on [0, L] N . 
The function iPl{x) := 4>l(Lx) thus satisfies —DAipL — L 2 w(x)ipL = L 2 Xi : li^l in [0, 1] N , with 
1-periodicity. From the Rayleigh formula we thus obtain 

. D /[o,lp IWP f [0 ,i]xVi> 2 
M,L — mm —r — t> -5 f jy'i 

therefore A^x < (since ifi = 1 G Hp er ), and Ai ; £ decreases with L. It implies that 62 increases 
with L. The relationship between Si and L is less clear since </>l = mine 4>L may not always be 
an increasing function of L. 

In order to lessen the boundary effects and to focus on fragmentation, we place ourselves in 
the periodic case. For our numerical computations, we assume that the environment is made 
of two components, favorable and unfavorable regions. This is expressed in the model (|2.ip 
through the coefficient n(x), which takes two values /j + or n~ , depending on the space variable 
x. We also assume that 

/i + > [i~ , v(x) = 1, h{x) = 1, and D = 1. 

Using a stochastic model for landscape generation |34j . we built 2000 samples of binary 
environments, on the two-dimensional period cell C = [0, l] 2 , with different degrees of frag- 
mentation. In all these environments, the favorable region, where fJ,(x) = occupies 20% of 
the period cell. The environmental fragmentation is defined as follows. We discretize the cell 
C into nc = 50 x 50 equal squares Cj. The lattice made of the cells Cj is equipped with a 
4-neighborhood system V(Q) (see Figure [2|), with toric conditions. On each cell Cj, we assume 
that the function [i takes either the value fi + or , while the number n+ = cardji, /j = fi + 
on Ci} is fixed to nc x ^ = 500. For each landscape sample ui = (/j(Cj))j = i v .. jnc , we set 
s(u) = I ^QeC Ylc-eV(Ci) ^-{^(Cj) = u(C{)}, the number of pairs of neighbors (Cj,Cj) such 
that /i takes the same value on Cj and Cj (!{•} is the indicator function). The number s(u>) is 
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Figure 2. The 4-neighborhood system: an element C{ of C and its four neighbors. 
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(a) s = 3400 (b) s = 3800 (c) s = 4200 (d) s = 4600 




10 1 



(e) s = 4800 (f ) s = 4900 

Figure 3. Some samples of the landscapes used for the computations of 5± and 
#2, with different values of the habitat aggregation index s. The black areas 
correspond to more favorable environment, where (J,(x) = 

directly linked to the environmental fragmentation: a landscape pattern is all the more aggre- 
gated as s(u) is high, and all the more fragmented as s(u)) is small (Figure [3]). Thus, we shall 
refer to s as the "habitat aggregation index." 

Remark 3.2. There exist several ways of obtaining hypothetical landscape distributions. The 
commonest are neutral landscape models, originally introduced by Gardner et al. |16j . They 
can include parameters which regulate the fragmentation [20] . We preferred to use a stochastic 
landscape model presented in [34], since it allows an exact control of the favorable and unfa- 
vorable surfaces and is therefore well adapted for analyzing the effects of fragmentation per 
se. This model is inspired from statistical physics. The number of pairs of similar neighbors 
s is controlled during the process of landscape generation. This quantity can be measured a 
posteriori on the landscape samples. Other measures of fragmentation could have been used, 
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1.3 



1.2 



1.1 



1 

3500 4000 4500 5000 

s: habitat aggregation index 

Figure 4. Solid lines: 8\j and $2j correspond respectively to the data sets 
{( s *> <5i)}i=i,.--,2000 and {(s l , ^2)}i=i,.--,2000) fitted with ninth degree polynomials. 
Dashed lines: 8\i Q is a lower prediction bound for new observations of 8\, and 
$2,up & n upper prediction bound for new observations of 82, with in both cases a 
certainty level of 99%. 

such as fractal dimension (see [M])- For a discussion on the different ways of measuring habitat 
fragmentation in real- world situations, the interested reader can refer to [14J. 

For our computations, we took fi + = 10 and fi~ = 0, and we computed the corresponding 
values of X\, 8\, and 5\ on each landscape sample ui' L of aggregation index s l , for i = 1, . . . , 2000. 
The eigenvalues X\ were computed with a finite elements method. We fitted the data sets 
{(s J , #i)}i=i,...,2000 and {(s l , ^2)}i=l,—,2000 using ninth degree polynomials (it is enough to assess 
whether the relations between s and 81,82 tend to be monotonic or not). The resulting fitted 
curves 8\j and 82 j are presented in Figure 01 Under the assumption of normally distributed 
values of 8% and 82 for fixed s values, we computed a lower prediction bound (<5i,; ) for new 
observation of 8\ and an upper prediction bound for 82 (82,up), with a level of certainty of 99%. 
Thus, given a configuration oj, with a fixed value of s, when 8 is smaller than 8\^ we take 
a 0.5% chance of b eing above <5i, while when 5 is larger than $2,up 

we take a 0.5% chance 

of being below 82- The small thickness of the intervals (8i t i Q , 82^) emphasizes the quality of 
the relationship between the habitat aggregation index s and the maximum sustainable yield 
8* E [<$i,<$2]- This also indicates that the criteria of Theorems 12.51 and 12.61 are close to being 
optimal, at least in some situations. 

Furthermore, as we can observe, the values of 81 and 82 tend to increase as s increases, 
and thus as the environment aggregates. Since 8* G [^1,^2], we deduce from the computations 
presented in Figure U] that 8* tends to increase with environmental aggregation. 

These tests were performed for particular values of /i + and [i~ . However, the thickness 
of the interval (81,82) can be determined for all values of fj, + ,fj,~ without further numerical 
computations, provided that fi + — pT = 10. Indeed, let us set B := fi + — . For a fixed value 
of B, let Ho(x) be a given L-periodic function in L°°(W N ) taking only the two values = B 
and Hq = 0. Let Ai 5 o be the first eigenvalue of the operator —V 2 — hqI on C, with L-periodicity 




16 



L. ROQUES, AND M.D. CHEKROUN 



conditions, 4>q the associated eigenfunction with minimal value (f>o, and 



81,0 ■= M ,' 7 \ o and 8 2 ,o 



(1 + 



We have the following proposition. 

Proposition 3.1. Assume that p(x) = fio(x) + wi/i /U~ > A^o- £e£ 5i and 82 be defined by 
([2~15]) . T/ien we force 5 2 - fa = (1 - ^) 2 (^2,o - «5i,o)- 

This result also indicates that the information on 8* is all the more precise as the growth rate 
function takes low values. However, the "relative thickness" of the interval (61,62), compared 
to 61, S2 ^ Sl , does not depend on fj,~, as can be easily seen. 

Proof of Proposition 13.11 The relation Ai[ / u(x)] = A^q — pT is a direct consequence of the 
uniqueness of the first eigenvalue Ai. We assume that pT > A^o, so that Ai[/i(x)] < 0. From 
the uniqueness of the eigenfunction (ft associated with Ai, <f> does not depend on \x~ . Therefore, 

6\ and 6 2 satisfy 8\ = ~7jq^rp ~ an d 82 = ^ Al '° 4 ^ - . The result immediately follows. □ 

4. A FEW COMMENTS ON THE PROPORTIONAL HARVESTING MODEL 
In this model, the population density u is governed by the equation 
(4.1) u t = DV 2 u + u(n(x) — u(x)u) — q(x)u, xefi, 

with L-periodicity of the functions p(x), u(x), and q(x) in the periodic case, and with Neumann 
or Dirichlet boundary conditions in the bounded case. Setting 

t(x) := n(x) - q(x), 

this model becomes equivalent to the SKT model (|1.3p . Hence, many properties of the solutions 
of this model are described in the existing literature. In particular the existence, nonexistence, 
and uniqueness results of Theorems 12.11 and 12.21 apply. The condition Xi[fj,(x) — q(x)] < is 
therefore necessary and sufficient for species persistence. Furthermore, the theoretical results of 
[8], [12], [33], [3l] on the effects of habitat arrangement on species persistence are also true for 
this model. 

For instance, when the function fi(x) is constant, with fi(x) = p\ > 0, and if the domain 0, 
is convex and symmetric with respect to each axis {x\ = 0}, . . . , {xn = 0}, the next result is a 
straightforward consequence of the paper [8]. 

Theorem 4.1. (i) In the periodic case, X\\pLi — q%(x)] < \\[pi — q(x)]. 

(ii) In the bounded Dirichlet case, \\\pi — ql(x)] < Ai[/ii — q(x)]. 

(iii) In the bounded Neumann case, if Q is a rectangle, Ai[/xi — g'l(ic)] < Ai[^i — q(x)]. 

Here q^. denotes the symmetric decreasing Steiner rearrangement of the function q with respect 
to the variable Xk, and q\ denotes the monotone rearrangement of q with respect to Xk (see [8j 
and [10] for the definition of these rearrangements). These rearrangements of a function q 
preserve not only its mean value, but also its distribution function. This means that if, for 
instance, q corresponds to a "patch" function taking the values q±, q2, and q% in some regions 
A\, A2, and A3, respectively, with A± + A2 + A3 = \C\, then the areas of the regions where the 
rearranged functions q* and take the values q\, q 2 , and q^ remain equal to A%, A 2 , and A3, 
respectively. 

Theorem 14.11 combined with Theorem 12.21 says that the spatially rearranged harvesting strate- 
gies are better for species survival. This result can be helpful from a resource management point 
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Figure 5. Examples of applications of Theorem l4.H ii-iii to reserves manage- 
ment. In panels (a) and (b), the boundary r of 0, is lethal (Dirichlet boundary 
conditions), (a) The initial effort function q{x) takes two values, q + > in the 
white area, and q~ = in the shadowed regions, which correspond to reserves, 
(b) Position of the reserves after a symmetric decreasing Steiner rearrangement 
along the Ai and A2 axes, successively. The rearranged configuration (b) always 
give more chances of species persistence. In panels (c) and (d), the boundary 
T is divided into two parts: r = Ti U IV Ti is represented with a solid line 
and can correspond to a coast, while T2 is represented with a dashed line and 
can correspond to a nonphysical limit that the species cannot cross (Neumann 
boundary conditions), (c) The effort function q(x) again takes two values, q + > 
in the white area, and q~ = in the reserves, (d) Position of the reserves af- 
ter monotone rearrangement along the horizontal and vertical axes, successively. 
The chances of persistence are better in the rearranged configuration (d). 
of view. Indeed, the authorities can rearrange the position of the harvested areas in order to 
improve the chances of population persistence. The result of Theorem 14.11 shows that, in the 
framework of these models, the creation of a large reserve gives persistence more chances than 
the creation of several small reserves, and is in accordance with the former results of |23j and [26] 
in the Dirichlet case. See Figure [5] for some illustrations in the bounded case with Dirichlet and 
Neumann boundary conditions. 
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5. Discussion 

We have proposed a model for the study of populations in heterogeneous environments, for 
populations submitted to an external negative forcing term. This forcing term could be regarded 
as a "quasi-constant-yield" harvesting, depending only on the population density u when u is 
below a certain small threshold e. The introduction of such a threshold e was necessary for 
ensuring the nonnegativity of the solutions of our model, and therefore its actuality. 

We carried out new mathematical results on the elliptic equation satisfied by the stationary 
states of the model, and on the associated parabolic equation. Both qualitative and quantitative 
results were obtained. 

From the qualitative point of view, we described the behavior of the model solutions in terms 
of the harvesting amplitude 5. Two main types of stationary solutions were found: the remnant 
solutions, always below a small threshold Eq an d therefore close to 0, and the significant solutions, 
always above this threshold, thus ensuring a time-constant yield. We discussed the maximum 
number of significant stationary solutions, which we found equal to 2, under a hypothesis of 
positivity of the second eigenvalue A2 of a linear operator. We further investigated the long-time 
behavior of the solution of our model, starting from a nonharvested population at equilibrium. 
We found a critical value 5* of the harvesting term amplitude, below which the population 
density tends over time to a significant stationary solution, and above which it converges to a 
stationary solution which is not significant. We also established quantitative formulae for some 
lower and upper bounds for 5*: 6\ and 62 , respectively. The threshold 82 has the additional 
property that, whenever the amplitude 5 is above 62, the population density decreases to a 
remnant stationary solution. 

The quantitative aspects of our study mainly consisted of discussing the effect of environ- 
mental fragmentation on these thresholds 5i and 62, and therefore on the interactions between 
environmental fragmentation and maximum sustainable yield. Namely, when computing the val- 
ues of 5\ and 82 on 2000 samples of stochastically obtained patchy environments, with different 
levels of fragmentation, we found an increasing relationship between these two coefficients and 
an environmental aggregation index s. This indicates that, for given areas of favorable and unfa- 
vorable regions, the harvesting quota that a species can sustain, while ensuring a time-constant 
yield, is higher when the favorable regions are aggregated. 

The reader may note that, in our model, the species mobility was not affected by the environ- 
mental heterogeneity. Such a dependence could be modeled by using a more general dispersion 
term, of the form V • (A(x)Vu), instead of DX7 2 u, where A(x) stands for the diffusion matrix (see 
[8], [36]). In fact, most of our results still work when the matrix A is of class C ,a (with a > 0) 
and uniformly elliptic, i.e., when there exists r > such that A{x) > tIn for all x £ £1. Indeed, 
Theorems 12. 1\ \2.2\ \2A\ 12.51 and 12.61 remain true under this more general assumption. However, 
the effects of environmental heterogeneity may differ, depending on the way A{x) and /J,(x) are 
correlated (see [U]). In the proportional harvesting case, the results of section 4 on the effects 
of the arrangements of the harvested regions may also not be valid with this dispersion term. 
However, in situations where A{x) takes low values (slow motion) when q(x) is low ("reserves"; 
see section 4), as underlined in [33], a simultaneous rearrangement of the functions A(x) and 
q(x) would lead to lower Ai values and therefore to higher chances of species survival. 
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